/*
 * lib.h
 *
 *  Created on: 2010.06.25
 *      Author: Tadas
 */

#ifndef LIB_H_
#define LIB_H_

#include <cmath>

#define GOLD 1.618034
#define GLIMIT 100
#define TINY 1e-20

#define INGOLD 1.0/GOLD
#define CIN (1.0 - INGOLD)

#define ITMAX 100


using namespace std;

double sign(double, double);
double findMax(double, double);
void shift2Points(double &, double &, double &);
void shift3Points(double &, double &, double &, double &);

void minbrak(double &ax, double &bx, double &cx, double &fa, double &fb, double &fc, double(*func)(double))
{
	double ulim, u, r, q, fu, dum;
	fa = func(ax);
	fb = func(bx);

	if (fb>fa)
	{
		shift3Points(dum, ax, bx, dum);
		shift3Points(dum, fb, fa, dum);
	}

	cx = bx+GOLD*(bx-ax);
	fc = func(cx);

	while (fb>fc)
	{
		r = (bx-ax)*(fb-fc);
		q = (bx-cx)*(fb-fa);;
		dum = findMax(fabs(q-r), TINY);
		u = bx-((bx-cx)*q-(bx-ax)*r)/2.0*sign(dum, q-r);
		ulim = bx+GLIMIT*(cx-bx);

		if ((bx-u)*(u-cx)>0.0)
		{
			fu = func(u);
			if (fu<fc)
			{
				ax = bx;
				bx = u;
				fa = fb;
				fb = fu;
				return; //tarp bx ir cx
			}
			else if (fu>fb)
			{
				cx = u;
				fc = fu;
				return; //tarp ax ir u
			}
			u = cx+GOLD*(cx-bx);
			fu = func(u);
		}
		else if ((cx-u)*(u-ulim) > 0.0)
		{
			fu = func(u);

			if (fu < fc)
			{
				dum = u+GOLD*(u-cx);
				shift3Points(bx, cx, u, dum);
				dum = func(u);
				shift3Points(fb, fc, fu, dum);
			}
		}
		else if ((u-ulim)*(ulim-cx) >= 0.0)
		{
			u = ulim;
			fu = func(u);
		}
		else
		{
			u = cx+GOLD*(cx-bx);
			fu = func(u);
		}
		shift3Points(ax, bx, cx, u);
		shift3Points(fa, fb, fc, fu);
	}
}

double findMax(double a, double b)
{
	if (a>b) return a;
	else return b;
}

double sign(double a, double b)
{
	if (b>=0) return a;
	else return -a;
}

void shift2Points(double &a, double &b, double &c)
{
	a = b;
	b = c;
}

void shift3Points(double &a, double &b, double &c, double &d)
{
	a = b;
	b = c;
	c = d;
}



double golden(double ax, double bx, double cx, double(*func)(double), double tol, double &xmin, int &n)
{
	double f1, f2, x0, x1, x2, x3, dum;
	x0 = ax;
	x3 = cx;

	if (fabs(cx-bx) > fabs(bx-ax))
	{
		x1 = bx;
		x2 = bx+CIN*(cx-bx);
	}
	else
	{
		x2 = bx;
		x1 = bx - CIN*(bx-ax);
	}

	f1 = func(x1);
	f2 = func(x2);

	while (fabs(x3-x0) > tol*(1+fabs(x1)+fabs(x2)))
	{
		n++;
		if (f2 < f1)
		{
			dum = INGOLD*x2+CIN*x3;
			shift3Points(x0,x1,x2,dum);
			dum = func(x2);
			shift2Points(f1, f2, dum);
		}
		else
		{
			dum = INGOLD*x1+CIN*x0;
			shift3Points(x3, x2, x1, dum);
			dum = func(x1);
			shift2Points(f2,f1,dum);
		}
	}
	if (f1 < f2)
	{
		xmin = x1;
		return f1;
	}
	else
	{
		xmin = x2;
		return f2;
	}
}

//---------------------------------------------------------------------------
double brent(double ax, double bx, double cx, double(*func)(double), double tol, double &xmin, int &n)
{
	int iter;
	double a, b, u, v, w, x, xm;
	double fu, fv, fw, fx;
	double e=0.0, d, etemp, p, q, r, toli, tol2;

	if (ax < cx)
	{
		a = ax;
		b = cx;
	}
	else
	{
		a = cx;
		b = ax;
	}

	x = w = v = bx;
	fx = fw = fv = func(bx);

	for (iter = 1; iter<ITMAX; iter++)
	{
		xm = 0.5*(a+b);
		tol2 = 2.0*(toli=tol*fabs(x)+tol);

		if (fabs(x-xm) <= (tol2-0.5*(b-a)))
		{
			xmin = x;
			n = iter;
			return func(x);
		}
		if (fabs(e) > toli)
		{
			r = (x-w)*(fx-fw);
			q = (x-v)*(fx-fw);
			p = (x-v)*q-(x-w)*r;
			q = 2.0*(q-r);

			if (q > 0) p = -p;

			q = fabs(q);
			etemp = e;
			e = d;

			if (fabs(q) >= fabs(0.5*q*etemp) || p <= q*(a-x) || p >= q*(b-x))
			{
				if (x >= xm) e = a-x;
				else e = b-x;
				d = CIN*e;
			}
			else
			{
				d = p/q;
				u = x+d;
				if ((u-a) < tol2 || (b-u)<tol2) d = sign(toli, xm-x);
			}
		}
		else
		{
			if (x > xm) e = a-x;
			else e = b-x;
			d = CIN*e;
		}

		if (fabs(d) >= toli) u = x+d;
		else u = x + sign(toli, d);
		fu = func(u);

		if (fu <= fx)
		{
			if (u >= x) a = x;
			else b = x;
			shift3Points(v, w, x, u);
			shift3Points(fv, fw, fx, fu);
		}
		else
		{
			if (u < x) a = u;
			else b = u;
			if ((fu <= fw) || (w==x))
			{
				v = w;
				fv = fw;
				w = u;
				fw = fu;
			}
			else if (fu <= fv || v==x || v==w)
			{
				v = u;
				fv = fu;
			}
		}
	}
	printf("Per daug zingsniu funkcijoje brent\n");
	xmin = 0.0;
	return 0.0;

}



#endif /* LIB_H_ */
